


readFile = function(rep, sim, nChr, nSnp, map, gvar, ntraits, training){
    fileName = paste("results/results", sim, rep, nChr, nSnp, map, gvar, training, ntraits, sep = ".")
    values = read.table(fileName)

    return(data.frame(sim=sim, nChr = nChr, nSnp = nSnp, map = map, gvar = gvar, ntraits = ntraits, training = training, impAcc = mean(values[,1])))

}

read_replicate = function(...) {
    base = readFile(rep = 1, ...)
    data = base$impAcc[1]
    count = 1
    for(rep in 2:10){
        try({

            tmp = readFile(rep, ...)
            data = data + tmp$impAcc[1]
            count = count + 1 

        })  
    }
    base$impAcc[1] = data/count
    return(base)

}

getValues = function(sim, chr, nSnp, map, gvar, ntraits, training){
    values = list()
    for(s in sim){
        for(ch in chr){
            for(ma in map){
                for(snp in nSnp){
                    for(g in gvar){
                        for(nt in ntraits){
                            for(tr in training) {
                                try({
                                    values[[paste(s,ch, ma, snp,g,nt,tr)]] = read_replicate(sim = s, nChr = ch, nSnp = snp, map = ma, gvar = g, ntraits = nt, training = tr)
                                })
                            }
                        }
                   }
                }
            }
        }
    }
    return(do.call("rbind", values))
}



panelColors = colorRampPalette(c("blue", "red"))(6)

makePlot = function(values, dep, x_label, ...){
    #Fix number of traits, change gvar.
    plot(c(), ylab = "", xlab = "", xaxt="n", ...)
    title(xlab = x_label, line = 2)
    title(ylab = "Imputation Accuracy", line = 2.5)

    x_vals = unique(values[[dep]])

    odds = which(1:length(x_vals)%%2 == 1)
    evens = which(1:length(x_vals)%%2 == 0)
    axis(1, at=x_vals, labels = rep("", length(x_vals)))
    axis(1, at=x_vals[odds], labels=x_vals[odds], tick = FALSE)
    axis(1, at=x_vals[evens], labels=x_vals[evens], tick = FALSE)


    nTraits = c(1, 5, 10, 25, 50, 100)
    for(i in 1:length(nTraits)){
        ntrait = nTraits[i]
        subData = values[values$ntraits == ntrait,]
        points(subData[[dep]], subData$impAcc, type = "b", col = panelColors[i])
    }
}
scale = .75
dev.new(width=12*scale, height = 6*scale)
par(mfrow = c(2, 3), cex = 1, las=1, bty="n", mar =  c(3, 4, 1, 1) + 0.1)

# Genetic Variance
values = getValues(sim = "base", chr = 10, nSnp = 500, map = 1, gvar=c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7), ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
makePlot(values, "gvar", x_label = "Heritability", xlim = c(0.1, .7), ylim = c(0, 1))

# Training size
values = getValues(sim = "base", chr = 10, nSnp = 500, map = 1, gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = c(100, 250, 500, 750, 1000))
makePlot(values, "training", x_label = "Training population size", xlim = c(100, 1000), ylim = c(0, 1))

# HD_markers
values = getValues(sim = "hd_markers", chr = 10, nSnp = c(100, 500, 1000), map = 1, gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
makePlot(values, "nSnp", x_label = "Number of markers", xlim = c(100, 1000), ylim = c(0, 1))

# Number of chromosomes
values = getValues(sim = "num_chr", chr = c(5, 10, 15, 20), nSnp = 500, map = 1, gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
makePlot(values, "nChr", x_label = "Number of chromosomes", xlim = c(5, 20), ylim = c(0, 1))

# Map
values = getValues(sim = "map_len", chr = 10, nSnp = 500, map = c(0.5, 1, 2, 4), gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
makePlot(values, "map", x_label = "Genetic map length (100cM)", xlim = c(0.5, 4), ylim = c(0, 1))

plot.new()
nTraits = c(1, 5, 10, 25, 50, 100)
legend("center", legend = nTraits, fill=panelColors, title = "Number of traits", ncol = 1, bty = "n")

dev.print(pdf, 'Figure1.pdf')




get_accuracy = function(chr = 10, nSnp = 500, map = 1, gvar = 0.5, ntraits = 100, training = 0){

    data_point = values[values$nChr == chr &
                        values$map == map &
                        values$gvar == gvar &
                        values$ntraits == ntraits &
                        values$training == training, "impAcc"]
    print(round(data_point, 3))
}

values = getValues(sim = "base", chr = 10, nSnp = 500, map = 1, gvar=c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7), ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
get_accuracy(gvar = 0.1, ntraits =1)
get_accuracy(gvar = 0.5, ntraits =1)
get_accuracy(gvar = 0.1, ntraits =100)
get_accuracy(gvar = 0.5, ntraits =100)


values = getValues(sim = "base", chr = 10, nSnp = 500, map = 1, gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = c(0, 100, 250, 500, 750, 1000))

get_accuracy(training = 100)
get_accuracy(training = 500)
get_accuracy(training = 0)



values = getValues(sim = "num_chr", chr = c(5, 10, 15, 20), nSnp = 500, map = 1, gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
get_accuracy(chr = 5)
get_accuracy(chr = 10)
get_accuracy(chr = 20)

# Map
values = getValues(sim = "map_len", chr = 10, nSnp = 500, map = c(0.5, 1, 2, 4), gvar=.5, ntraits=c(1, 5, 10, 25, 50, 100), training = 0)
get_accuracy(chr = 1)
get_accuracy(chr = 2)
get_accuracy(chr = 4)

